% GAUTIER LE BIHAN - 2020
% Replication files for "Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector Menu-Cost Model?" Review of Economics and Statistics
%
% Table 5
clear;
tic


addpath('..\..\Utilities')  

load actual_moments_k

max_sectors=40;
load actual_moments_k
weight_j=actual_moments_k(:,12)/100;

load ..\..\Estimation_param\MS_produits_CPlus\param
param=param;
%*******************************************STAT SECTEURS PARAM
coicop1=actual_moments_k(:,1);
weight=actual_moments_k(:,12);
sum(weight);
sample_all=[coicop1 weight param];

moments_mean=(sample_all(:,9:end)'*weight)./sum(weight);
size(sample_all);

%/************* MEDIAN ALL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%TOTAL
p0=abs(tanh(sample_all(:,3)));
mu_c=abs(sample_all(:,4));
sig_a=exp(sample_all(:,5));
rho_a=ones(227,1)*0.6946;

param_all=[p0 mu_c sig_a rho_a];
weight=sample_all(:,2);
v=size(param_all,1);
weight_all=weight;

wmean_all=(param_all'*weight_all)./sum(weight_all);

%%%%%%%%%%%%%%%%%MEDIAN ALL EXCEPT GASOLINE
cond=(sample_all(:,9)>=0.5);
sample_e=sample_all((cond==1),:);
%%%%%%%%%%%%%%%%%%%%%%TOTAL
p0=abs(tanh(sample_e(:,3)));
mu_c=abs(sample_e(:,4));
sig_a=exp(sample_e(:,5));
rho_a=ones(size(sample_e,1),1)*0.6946;

param_e=[p0 mu_c sig_a rho_a];
weight=sample_e(:,2);
v=size(param_e,1);
weight_e=weight;


wmean_e=(param_e'*weight_e)./sum(weight_e);

%condition  sur dist de l'estimation 10:OK 221 secteurs 5 : 217 secteurs
cond=(sample_all(:,9)<10000);
sample=sample_all((cond==1),:);


size_p=size(sample,1)
coicop1=sample(:,1);

moments_mean=(sample(:,9:end)'*sample(:,2))./sum(sample(:,2));

moments=sample(:,9:end);
weight=sample(:,2);

%sample_clean=test((test(:,13)<1),:);
%save sample_clean sample_clean
%%%%%%%%%%%%%%%%%%%%%%TOTAL
p0=abs(tanh(sample(:,3)));
mu_c=abs(sample(:,4));
sig_a=exp(sample(:,5));
rho_a=ones(size(sample,1),1)*0.6946;

param=[p0 mu_c sig_a rho_a];
weight=sample(:,2);
v=size(param,1);
weight_tot=weight;



wmean=(param'*weight_tot)./sum(weight_tot);

for jj=2:max_sectors;
    number=jj;
%figure;hist(param_sect(:,3));

percentiles=ones(number-1,5);
    for i=1:size(param,2)
        [sortx,order] = sort(param(:,i));
        sortw = weight(order);
        
        for n=1:number-1
         midpoint = n*sum(sortw)/number;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              percentiles(n,i) = mean(sortx([j j+1]));
            else
                
                percentiles(n,i) = sortx(j+1);
            end
     percentiles(n,end)= csumw(end);
        end
    end  
    all_per(jj-1).all_per=percentiles;
end 

for zz=1:1:3
%freq=moments(:,1);
choix_param=zz;
sig=param(:,choix_param);
quantile=zeros(size_p,max_sectors-1);

for jj=1:max_sectors-1

per_sig=all_per(jj).all_per;
per_sig=per_sig(:,choix_param);
var_dum=zeros(size_p,jj);
    for n=1:jj
    var_dum(:,n)=(sig>=per_sig(n));
   
    end
    if jj==1
        quantile(:,jj)=var_dum;
    else
     quantile(:,jj)=sum(var_dum')';
    end
end


for kk=1:max_sectors-1

quantile_k=quantile(:,kk);

     wmean_quant=zeros(max(quantile_k)+1,size(param,2));
     weight_quant=zeros(max(quantile_k)+1,1);
     for k=1:max(quantile_k)+1
    
     param_sect=param((quantile_k==k-1), :);
     v=size(param_sect,1);
     weight_sect=weight((quantile_k==k-1),:);
    
     weight_quant(k,:)= sum(weight_sect);
     wmean_quant(k,:)= (param_sect'*weight_sect)/weight_quant(k,:);
     end  

wmean_quant=[wmean_quant];
weight_quant=[weight_quant];

 param_ms(kk+1).param_ms=wmean_quant;
  weight_ms(kk+1).weight_ms=weight_quant;
end;

 param_ms(1).param_ms=[wmean_all'];
  weight_ms(1).weight_ms=57.0181;
  

s=40;
%stat_mean=zeros(s,7);
%for jj=[1 2 3 4 5 6 7 9 11 13 15 17 19];
 % for jj=[3 4 5 6 7 9 11 13 15 17 19];
 stat_sect=zeros(23,120);
  stat_aggall=zeros(23,7);
 IRF_all=zeros(201,23);
   for tt=5:4:17;
   
    param_med=param_ms(tt).param_ms;
    param_med=param_med((param_med(:,1)~=0),:)
    
    weight_sect=weight_ms(tt).weight_ms;
    weight_sect=weight_sect((weight_sect(:,1)~=0),:)
    weight_sect
    p0=param_med(:,1)';
    mu_c=param_med(:,2)';
    sig_eps_a=param_med(:,3)';
    rho_a=param_med(:,4)';   
    if tt==1
        w=1;
    elseif tt==2
        w=[weight_sect sum(weight_e)];
    else
     w=[weight_sect'];
     end
   
    sm_m=0;
    param_e=[p0 mu_c sig_eps_a rho_a];

    [avmoments, moments, VC, IRF]=geNCalvoPlus_2Agg(sm_m, p0, mu_c, sig_eps_a, rho_a, w);
 NSect = size(w,2);
 moments
  
 
    stat_sect1(tt,1:NSect*6)=[ moments]
    save stat_sect1 stat_sect1;

   
   end

    temp=[ stat_sect1(5,:);stat_sect1(9,:);stat_sect1(13,:);stat_sect1(17,:) ]
  st_sector1(zz).st_sector1= temp;
  save st_sector1 st_sector1;
   

end;
 
 
 
 
 